////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//  -------------------- Ebrahim Foulaadvand & Somayyeh Belbasi, 30 July 2009 --------------------------- //
//                                                                                                        //
//  The programme "ForcedCoupledOscillators" evaluates the temporal evolution of a system of forced       //
//  coupled harmonic oscillator which are connected to each other by springs. It also gives you the       //
// amplitude frequency dependence of oscillators. The fix boundary condition is used. The second order    //
// Runge-Kutta and Euler-Cromer algorithms are implemented for solving the motion equations.              //
//                                                                                                        //
//                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>

using namespace std;


main()
{

randomize();

const N=3, T=30000; // N= number of oscillators, T= number of simulation timesteps.

double delt=0.01,F0=0.5,w=1,delw=0.05,E; // delt=timestep, F0=driving force amplitude.

//w=Sinusoidal driving force frequency, delw=increment for w.

int i,j,s,S=10,t;

vector< double > m(N+2,0),k(N+3,0),D(N+2,0),xEC(N+2,0),vEC(N+2,0),x(N+2,0),v(N+2,0),k1x(N+2,0),k2x(N+2,0),
                 k1v(N+2,0),k2v(N+2,0);

// Arrays "m" and "k" store the oscillators masses and the spring constants. Arrays "x" and "v" store
// oscillators positions and velocities computed by RK2 algorithm. Arrays "xEC" and "vEC" store
// oscillators positions and velocities computed by EC (Euler-Cromer) algorithm. We use fixed boundary
// condition i.e.; end masses are only connected to one adjacent mass. The first mass is connected to
// the driving force.


ofstream file1 ("xEC 2.plt"); //output file for the poistion time series of particle 2 (EC algorithm).
ofstream file2 ("xRK 2.plt"); //output file for the poistion time series of particle 2 (RK2 algorithm).
ofstream file3 ("D1(w).plt"); //output file for the frequency dependence of particle 1.
ofstream file4 ("D2(w).plt"); //output file for the frequency dependence of particle 2.
ofstream file5 ("D3(w).plt"); //output file for the frequency dependence of particle 3.


for (int i=1;i<=N;i++){ m[i]=1; k[i]=1; }

k[N+1]=1;



//----------- Euler-Cromer algorithm --------------------



for (int s=1;s<=S;s++){ w=s*delw; cout<<"s= "<<s<<endl;

xEC[0]=0; vEC[0]=0; xEC[N+1]=0; vEC[N+1]=0;

for (int i=1;i<=N;i++){ xEC[i]=0; vEC[i]=0; D[i]=0; }


for (int t=1;t<=T;t++) {

vEC[1]=vEC[1]- delt*(k[1]/m[1])*xEC[1] - delt*(k[2]/m[1])*(xEC[1]-xEC[2]) + delt*F0*cos(w*t*delt);

for (i=2;i<=N;i++) {

vEC[i]=vEC[i]- delt*(k[i]/m[i])*(xEC[i]-xEC[i-1])-delt*(k[i+1]/m[i])*(xEC[i]-xEC[i+1]);

} // i


for (i=1;i<=N;i++) { xEC[i]=xEC[i]+ delt*vEC[i];  if ( fabs(xEC[i]>D[i]) ) D[i]=fabs(xEC[i]); }

file1<<t*delt<<" "<<xEC[2]<<endl; // storing in the data file


}

//----------- 2nd Order Runge-Kutta algorithm --------------------

x[0]=0; v[0]=0; x[N+1]=0; v[N+1]=0;

for (int i=1;i<=N;i++){ x[i]=0; v[i]=0; D[i]=0; }


for (int t=1;t<=T;t++) { //cout<<"t= "<<t<<endl;getch();


k1v[1]=-delt*(k[1]/m[1])*x[1] - delt*(k[2]/m[1])*(x[1]-x[2]) + delt*F0*cos(w*t*delt);
k1x[1]=delt*v[1];


k2v[1]=-delt*(k[1]/m[1])*(x[1]+0.5*k1x[1]) - delt*(k[2]/m[1])*(x[1]+0.5*k1x[1]-x[2]-0.5*k1x[2]) + delt*F0*cos(w*(t+0.5)*delt);
k2x[1]=delt*(v[1] +0.5*k1v[1]);


for (i=2;i<=N;i++) {

k1v[i]=-delt*(k[i]/m[i])*(x[i]-x[i-1])-delt*(k[i+1]/m[i])*(x[i]-x[i+1]);
k1x[i]=delt*v[i];


k2v[i]=-delt*(k[i]/m[i])*(x[i]+0.5*k1x[i]-x[i-1]-0.5*k1x[i-1])-delt*(k[i+1]/m[i])*(x[i]+0.5*k1x[i]-x[i+1]-0.5*k1x[i+1]);
k2x[i]=delt*(v[i] +0.5*k1v[i]);

} // i

for (i=1;i<=N;i++) {

v[i]=v[i] + k2v[i];
x[i]=x[i] + k2x[i];

if ( fabs(x[i]>D[i]) ) D[i]=fabs(x[i]);

}

file2<<t*delt<<" "<<x[2]<<endl; // storing in the data file


} // t

file3<<delw*s<<" "<<D[1]<<endl; // storing in the data file
file4<<delw*s<<" "<<D[2]<<endl; // storing in the data file
file5<<delw*s<<" "<<D[3]<<endl; // storing in the data file


//cout<<"D[1]= "<<D[1]<<endl;//getch();

} // end s


cout<<"end."<<endl;//getch();


return 0;

} //main